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In the numerical renormalization group (NRG) calculations of spectral functions of quantum impurity mod- 
els, the results are always affected by discretization and truncation errors. The discretization errors can be 
alleviated by averaging over different discretization meshes ("z-averaging"), but since each partial calculation 
is performed for a finite discrete system, there are always some residual discretization and finite-size errors. 
The truncation errors affect the energies of the states and result in the displacement of the delta peak spectral 
contributions from their correct positions. The two types of errors are interrelated: for coarser discretization, 
the discretization errors increase, but the truncation errors decrease since the separation of energy scales is en- 
hanced. In this work, it is shown that by calculating a series of spectral functions for a range of the total number 
of states kept in the NRG truncation, it is possible to estimate the errors and determine the error-bars for spectral 
functions, which is important when making accurate comparison to the results obtained by other methods and 
for determining the errors in the extracted quantities (such as peak positions, heights, and widths). The closely 
related problem of spectral broadening is also discussed: it is shown that the overbroadening contorts the results 
without, suiprisingly, reducing the variance of the curves. It is thus important to determine the results in the 
limit of zero broadening. The method is applied to determine the error bounds for the Kondo peak splitting in 
external magnetic field. For moderately large fields, the results are consistent with the Bethe Ansatz study by 
Moore and Wen. We also discuss the regime of large U /T ratio. It is shown that in the high-field limit, a spectral 
step is observed in the spectrum precisely at the Zeeman frequency until the field becomes so large that the step 
merges with the atomic spectral peak. 

PACS numbers: 72.10.Fk, 72.15.Qm 



I. INTRODUCTION 



Quantum impurity physics^ 2 - is an active area of research, 
which is particularly important for the problems of magnetic 
properties of confined electrons (quantum dots^i, magnetic 
impurity atoms on surfaces^), but also for strongly corre- 
lated electron systems due to the mapping of bulk correlated 
models to self-consistent single-impurity models within the 
dynamical mean-field theory^. For strong electron-electron 
interactions, the quantum impurity models are notoriously 
difficult to solve and they generally require the application 
of non-perturbative techniques. One such technique is the 
numerical renormalization group (NRG) 2 ' 9 ' 10 , which consists 
of discretizing the continuum of the conduction-band states, 
transforming the problem to the form of a semi-infinite tight- 
binding chain, and numerically diagonalizing the resulting 
discrete Hamiltonian in an iterative way. The discretiza- 
tion is performed by splitting the energy band into inter- 
vals of widths that decrease as a geometric series (oc A~ n , 
where A > 1 is called the discretization parameter) as the 
Fermi level is approached 9 . This particular choice of the dis- 
cretization scheme is adapted to the behavior of the Kondo 
model, where each energy scale makes a comparable con- 
tribution to the renormalization group flow of the exchange 
coupling constant in the scaling regime 9 ^. The NRG was 
first applied to calculate the thermodynamic properties of im- 
purity problem a 9 ' 10,12 ^ , and was later extended to dynami- 
cal properties^—. Further important improvements were the 
development of the density-matrix approach to spectral func- 
tion calculation 2 ^, the self-energy trick^, and the introduction 



of the complete-Fock-space basis* 2 -— which solved the over- 
counting problem. 

Since the calculation is performed for a discretized prob- 
lem, one expects significant systematic discretization errors. 
They appear, for example, in the form of oscillations in 
the calculated spectral functions with frequency log A and 
its harmonics (in logarithmic frequency space). These os- 
cillations can be reduced by performing the so-called z- 
averaging, wherein one performs the same NRG calculation 
for several interleaved discretization meshes and averages the 
results 2 ^—. By averaging over two meshes, one cancels out 
log A oscillations, by averaging over four meshes one cancels 
out 2 log A oscillators, etc., thus the z-averaging is best per- 
formed for N z meshes where N z = 2™. Using an improved 
discretization scheme 2 -, the cancellation of oscillations is re- 
markable even in the case of strong hybridization of the impu- 
rity with the conduction band states and for large values of the 
discretization parameter A. Nevertheless, the spectral func- 
tions calculated using the NRG are always affected by the 
discretization and the finite-size errors to some degree, even 
when all technical refinements are used 2 ^. 

Another source of systematic errors in the NRG is the trun- 
cation. Since the Fock space grows exponentially with the 
chain length (by a factor of 4 for a spinfull single-channel 
impurity problem), the set of states kept after each step is 
truncated to some finite number N. This is clearly an ap- 
proximation, which was, however, shown to lead to highly ac- 
curate results 9 . The approach works because of the "energy- 
scale-separation" property of quantum impurity problems: the 
matrix elements coupling high-energy and low-energy excita- 
tions are small and controlled by A, large A leading to stronger 
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decoupling. 

Finally, for particle-hole asymmetric baths there is a further 
source of error in the NRG (the "mass-flow effect") due to the 
iterative algorithm used in the NRG for integrating out the 
impurity bath degrees of freedom, since the impurity parame- 
ter shift due to the real part of the bath propagator at a given 
NRG step only includes the contribution from the chain sites 
already included in the calculation, while the the contribution 
of the remaining half-infinite chain is missing 2 ^. The mass- 
flow effect is particularly severe for bosonic baths, while for 
fermionic baths it was found to have little effect^ 9 .. 

While the presence of the systematic errors in the NRG cal- 
culations is common knowledg e) 2 ' 9 ' 29 , few systematic studies 
have actually appeared in the literature and the dependence 
of the errors on the calculation parameters is still not widely 
known. The purpose of this work is to analyze the discretiza- 
tion and truncation errors in the spectral functions obtained in 
the most sophisticated calculations using the complete Fock- 
space (CFS) approach 24,25 with the self-energy trick 21 and the 
improved discretization mesh with averaging over many val- 
ues of the z -parameter 28 . It will be shown that with the in- 
creasing number N of states kept, the spectral function ob- 
tained in the CFS approach exhibits variations which arise 
due to the discrete nature of the NRG Hamiltonian and the 
particular way of collecting the spectral information in the 
CFS technique 2 ^, which is susceptible to truncation errors. By 
calculating the statistical properties of these variations, one 
can obtain useful information which quantifies the unavoid- 
able discretization and truncation errors in the NRG. The er- 
ror estimates thus obtained are actually lower error bounds, 
since it is conceivable that in addition to the errors which 
lead to variations as a function of the NRG calculation pa- 
rameters (N, A, and the number of points in the z-averaging) 
there are other systematic errors which do not average out, 
thus the NRG spectra may deviate more from the true spec- 
tra than the proposed error estimates indicate, but there is no 
way to detect such effects within the NRG itself. Neverthe- 
less, the knowledge about the amplitude of the oscillatory dis- 
cretization and truncation errors, especially in relation with 
the spectral broadening problem (discussed below), is impor- 
tant to make the best possible use of the method. 

H. MODEL AND DETAILS OF THE NUMERICAL 
TECHNIQUE 

We study the single-impurity Anderson model 1 ' 10 H = 
+ H c + H hath with 

-"imp — 

cr 

while H c and ffbath are the coupling and the conduction band 
parts. The electron repulsion is U = 0.1, the on-site en- 
ergy is Ed = —U/2, the hybridization is T = 0.008, and 
the external magnetic field is b = gj-LsB = 3.273 10~ 4 . 
All parameters are expressed in units of the half-bandwidth 
D = 1. Unless otherwise noted, the discretization parame- 
ter is A = 2, N z = 12 discretization meshes are used, and 



the broadening parameter is a ~ 0.075 using the broadening 
kernel proposed in Ref. (After rescaling by D = 10 and 
adopting a different convention for expressing the hybridiza- 
tion strength as TlD = 0.16, these are the same parameters 
as used in Ref. |3j]; in a later section, we will, as an exam- 
ple, provide an approximation for the Kondo-peak splitting 
together with an error estimate for this parameter set.) The 
spectra are broadened and z-averaged before the self-energy 
trick is applied. All calculations have been performed for 
zero temperature; in this case the CFS method 24 and the full 
density-matrix method 2 - become fully equivalent. In the NRG 
calculation, the same maximum number of states N is kept 
in the truncation for all values of z; this is necessary for a 
meaningful z-averaging in the CFS approach. (If, instead, the 
spectral functions are calculated using the alternative patching 
approac h 19 ' 32 , it is advantageous to use the truncation with a 
fixed energy cutoff; see also Ref.[28l where the systematic er- 
rors in the patching approach are comprehensively analyzed). 
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Figure 1: Spectral function A^(uj) of the Anderson impurity model 
calculated for a range of N, the number of states kept in the trunca- 
tion after each NRG iteration step. N ranges from 1800 to 3600 in 
steps of 100. The Kondo temperature (defined as in Refs. I MToh is 
T K = 6.9 10~ 5 , thus B/T K = 4.7. 

In Fig. Q] we plot the spectral functions obtained in the cal- 
culations for a variable number of states kept, ranging from 
N = 1800 (the value used in Ref. [0) up to N = 3600. It 
has to be noted that even the lower limit N = 1800 is suf- 
ficient to obtain essentially fully converged NRG results for 
static quantities (ground state energy, thermodynamic func- 
tions); the plot in Fig. [2] shows that at N = 1800 the error 
in the calculated ground state energy is below 10~ 6 , which is 
already a highly precise result. As evidenced in Fig. [1] the 
spectral functions, however, do not converge to some limiting 
curve as N is increased. This is due to the way the spectral 
function is calculated in the CFS approac h 24,25 : the delta-peak 
contributions appear at frequencies uj which are a difference 
of the energy of a kept state Ef. and the energy of a discarded 
state Ed, that is, uj = Ed — E^. While the kept states are in 
the part of the on-shell excitation spectrum which is expected 
to be rather accurate, thus Ek is precise, the truncated states 
come from the top of the spectrum which is more significantly 
affected by the accumulated truncation errors from the previ- 
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ous NRG steps, thus the energies Ed of these states are known 
to a lesser precision. In other words, in the CFS approach the 
normalization of the spectral function is guaranteed to be ex- 
actly 1, but the higher moments are not exact, i.e., only the 
0-th spectral sum rule is fulfilled to numeric precisio n 24,28 . In- 
creasing N does not help in this respect, since this merely 
implies that a given excitation will contribute at a later NRG 
step, thus Ed may accumulate even more truncation error. For 
this reason, changing N will change the resulting spectrum. 
Since the states in the NRG are clustered, the changes can be 
relatively abrupt. (We note in passing that the spectral func- 
tions calculated using the patching approach converge as N 
is increased, because one extracts the spectral information al- 
ways from the same energy interval of the on-shell spectra, 
thus the effect of the truncation errors of the discarded states 
is tiny. Alas, the patching approach suffers from other defi- 
ciencies; in particular, there is a free parameter which needs 
to be tuned for each particular application, which limits the 
reliability of the approach.) 
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Figure 2: (Color online) Ground-state energy _Egs as a function of 
the number of states kept, N. Eas is computed as described in 
RefJH- EqI is the ground-state energy of the conduction band with- 
out the impurity, thus what is plotted is the impurity binding energy. 
The variation of the binding energy (in the large- N limit) with A is a 
direct consequence of the discretization errors. These results indicate 
the degree of convergence of the NRG calculation as a function of iV 
and A. 



IE. ESTIMATION OF THE SPECTRAL-FUNCTION 
VARIANCE 

Fig. [3] shows the spectral function A(ui) obtained by aver- 
aging over the curves shown in Fig. Q] together with the con- 
fidence region which represents the estimated range of values 
over which the spectral function fluctuates as a function of the 
truncation cutoff N, The confidence region is determined by 
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Figure 3: (Color online) Upper panel: Spectral function calculated by 
averaging the CFS spectra over all values of N (thick black line) with 
a confidence interval determined by the standard deviation of the 
data at each uj (thinner lines/red online). For comparison, the spec- 
tral function calculated using the patching approach is also shown 
(dashed line/green online). The gray lines in the background are the 
individual CFS spectra from Fig.fTJ Lower panel: standard deviation 
as a function of frequency. 



calculating the standard deviation a(uj) at each frequency w, 
the upper and lower boundaries of the region are then taken to 
be A(uj) ± <t{oj). For comparison, the spectral function cal- 
culated using the traditional patching approach is also shown; 
it is evident that this curve lies within or near the confidence 
region of the spectral function calculated using the CFS ap- 
proach. In this sense, the two approaches appear equivalent in 
this case, however the patching approach by itself provides 
no means for estimating the reliability (i.e., the confidence 
region) of the result. The lower panel in Fig. [3] shows the 
standard deviation, cr(<j). It attains its highest values near the 
atomic peaks at uj ~ e<j, £rf + U and it tends to decrease at 
lower frequencies, as expected, although it has a small local 
maximum at the Kondo resonance where the results are again 
more scattered. 

From the results in Fig. [3] we extract the position of the 
Kondo resonance as the frequency of the maximum of the av- 
eraged spectral function, w max = 2.50 10 -4 . Let us consider 
now the estimation of the error committed due to the vari- 
ance of the spectral functions. As a simple (albeit pessimistic) 
approximation we may consider that the true maximum is lo- 
cated anywhere within the triangle plotted in Fig. |4] The av- 
erage width of the triangle is (wg — tiu)/2, thus the error 
estimate can be defined as 



From the calculations for the discretization A = 2 and the 
broadening a = 0.075 we thus conclude that the Kondo reso- 
nance is shifted by the external magnetic field to 

uj k w 2.50 10~ 4 ± 1.2 10~ 5 w 2.50 10" 4 (1 ± 0.05) , (3) 
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IV. SPECTRAL FUNCTION BROADENING 



w K /i)K 0.76(1 ±0.05), 



(4) 



where b = g^sB is the magnetic field in energy units (Zee- 
man energy). In other words, using the NRG calculations at 
A = 2, the position of the Kondo peak can be determined at 
best with 5% accuracy. This is the reason why it is so difficult 
to reliably study the dependence between the Kondo peak po- 
sition and the external magnetic field, which should behave as 
ujf = (2/3)6 at low fields (b < k B T K j££. and as u K = b 
at high fields^— (b > k B T K , but B still small compared to 
the atomic parameters of the model, otherwise non-universal 
features are observe d 31,41 ). 




Figure 4: (Color online) Estimation of the true position of the peak 
maximum; see the text for details. The data is taken from the left 
inset in Fig. [3] 

To gain more insight into the origin of the variation of the 
spectral function with N, in Fig. [3] we plot the raw binned 
spectral data for a fine-grained set of 256 different z values. A 
comparison of the results for N = 2000 and N = 3000 shows 
that while the general aspects do not vary with N, the details 
do. The differences in these details lead to the variation of the 
spectral function with increasing N. 
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Figure 5: Raw binned spectral data in the energy range of the Kondo 
peak for a fine-grained range of z values (N, — 256). These are raw 
spectral weights of the delta peaks forming the spectral functions 
which have been binned (logarithmic mesh, 600 bins per decade; the 
intervals of the binning grid are much narrower than the broadening 
kernels used to post-processes these raw results). Left and middle 
panels show results for N — 2000 and N = 3000, while the right 
panel shows the difference between the two. The grayscale ranges 
are approximately equal in the three plots. In the plots z ranges from 
to 1 (bottom to top) and u> ranges from 0.0002 to 0.0003 (across 
the Kondo peak). The model parameters are as in Fig.Q] 



We now discuss the role of the spectral function broadening 
parameter a which is used to obtain smooth spectral curves 
from the spectral information in the form of weighted delta 
peaks. The broadening kernel used is 



P(lu,lo') = K > exp 
\/7ra|w| 



\og\u/(jj'\ a 
a 4 



(5) 



i.e., the broadening kernel proposed in Ref . l25l with 7 = a/4. 
The peak of this kernel function is located at 



u = e- a /V. 



(6) 



The weight is distributed asymmetrically with respect to u>' 
with 1/2[1 + erf (a/4)] of the weight in the interval \u)\ > \u'\ 
and 1/2[1 — erf(a/4)] in the interval |w| < |u/|, where erf 
is the error function. Furthermore, the tail of the broadening 
function reaches to relatively high frequencies if a is large, 
thus the presence of high-energy spectral features might lead 
to spurious shifts of the peak positions in the low-energy part 
of the spectrum, and vice versa. 

Naively, one could expect that increasing a will reduce the 
variance of the spectral curves. This is only partially true: the 
irregularities indeed smooth out, however the overall spectral 
weight distribution is still found to be fluctuating with N. Sur- 
prisingly, it is found that that Scjk, the estimated (pessimistic) 
error in the Kondo peak position, remains roughly constant 
with a, see Fig. [6] for A = 2. It is thus important to use a 
small broadening parameter to avoid systematic overbroaden- 
ing error a 31 ' 42 . Clearly, a suitably small a is such that the over- 
broadening error is smaller than the intrinsic error due to the 
discretization (for instance, no larger than a = 0.2 or a = 0.3 
in our example, see Fig.|6]l. The broadened spectral functions 
depend on the broadening kernel used; only in the limit of very 
small broadening widths do the different broadening kernels 
(Gaussian, log-Gaussian, modified log-Gaussian) all become 
equivalent. For moderate and large broadening, it was found 
that the modified log-Gaussian kernel works best, see Fig. 19 
in Ref.il 

It may be noted that the variance does not depend on the 
choice of the discretization scheme, because for the chosen 
model parameters and A = 2, the discretization artifacts are 
small in any discretization scheme. 

Alternatively, one can determine the error in the extracted 
Kondo peak position by locating the maximum for each spec- 
tral function calculated at fixed N, and calculating the stan- 
dard deviation. The results of such a calculation are shown 
in Fig. |7] The error bars are significantly smaller using this 
definition and they decrease with increasing a. We can de- 
scribe this procedure as "optimistic", since it is found that the 
spectral functions with different N in fact have their maxima 
at positions which fluctuate less that the curve does overall, 
thus the error bars are significantly smaller than in the "pes- 
simistic" case presented in Fig. [6] The true systematic er- 
ror probably lies somewhere between these two extreme er- 
ror estimates, thus a suitable value of a for performing calcu- 
lations aiming towards high-precision results is smaller than 
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than expected from the shift described by Eq. ©. A good fit 
to the results for A = 2 in Figs.[6]and[7]is 

u K (a) = WA -(0)e- Q2/& (7) 

with b w 2.25. Thus the functional form of the shift is similar 
to that of Eq. ©, however the numerical factor in the argu- 
ment of the exponential function is 2.25 rather than 4, i.e., the 
shift is even larger than expected. Test calculations for a sin- 
gle Lorentzian peak indeed show that spectral peaks with a 
finite width are broadened into a wider peak whose maximum 
is displaced by a factor of e~" / b with b smaller than 4. This 
is another reason for reducing a as far as possible. 



Figure 6: (Color online) The Kondo peak position (with "pes- 
simistic" error bars) as a function of the broadening parameter. The 
(average) peak position is extracted by determining the maximum 
of the average of the spectral functions for different N; the maxi- 
mum is located using the conjugate gradient method. The smallest 
a, for which the Kondo peak maximum can still be extracted, is de- 
termined by the number N z of different values of the twist parameter 
(here N z = 12, and a m in ~ 0.05). For completeness, the results 
for larger broadening parameters A = 3 and A = 4 are also shown; 
smaller error-bars result from the fact that for coarser bro adening the 
truncation errors decrease (while the discretization errors increase). 
Thus larger values of A actually lead to less scatter in the calculated 
spectral functions as N is varied. 



suggested by the error bars in Fig. |6] i.e., 
rameter should probably be closer to a 
a = 0.2 or 0.3, as suggested above. 
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Figure 7: (Color online) The Kondo peak position (with "optimistic" 
error bars) as a function of the broadening parameter. The peak po- 
sition is extracted as an average of the maxima of the spectral func- 
tions for different N; each maximum is located using the conjugate 
gradient method. The error bars are determined by calculating the 
standard deviation of the extracted values for ujk- The error bars 
increase significantly for a < 0.1 because at given number of the 
twist parameters, N z = 12, the oscillations in the spectral functions 
become severe for a < l/y/AN z . 

It is interesting to note that the effect of finite a is larger 



V. KONDO RESONANCE SPLITTING 

As an application, we now consider the problem of the 
Kondo resonance splitting in an external magnetic field. We 
plot the splitting ratio ujk / (9HbB) as a function of the ra- 
tio g/iBB/ksTK in Fig. [8] (Note that uik is defined as the 
shift of the Kondo peak position in the spin-projected spectral 
function A a (u>) and that the peak-to-peak distance in the spin- 
averaged spectral function is not exactly 2u>k-) Both proce- 
dures for extracting the average value and the error bars have 
been performed. The average values (circles and squares in 
the figure) agree within the error bars of the "optimistic pro- 
cedure" for all the results shown; deviation becomes larger 
for smaller magnetic fields. For large fields, the ratio in- 
creases in a rather slow (logarithmic) way, thus one expects 
non-universal features to appear before the universal high- 
field asymptotic behavior is reached. For small fields, the 
ratio goes toward a value of 0.7, which is close to the ex- 
pected low-field limit valu o 34 ' 35 of 2/3. [The low-field asymp- 
totic limit has been reported to be confirmed in a calculation 
where broadening is performed using Lorentzian peaks with 
constant width at very low energy scales 3 ^.] The "pessimistic" 
error bars grow larger for small fields, and the "optimistic" 
average values start to deviate from the "pessimistic" ones. 
This is expected, since the Kondo peak displacement becomes 
smaller than the Kondo peak width, thus the relative errors 
(i.e., the error of the ratio ujk I 9^bB) grow with decreasing 
B because the absolute error Sujk approximately saturates for 
g^ B B < k B T K . 

The results in Fig. [8] are in good agreement with the split- 
ting in the Kondo model as determined from the spinon den- 
sity of states in the Bethe Ansatz (BA) solution 3 -^ see the full 
line with crosses in Fig. [8] At large fields, the deviation be- 
yond the error bars is most likely due to the differences be- 
tween Anderson and Kondo models at high fields which leads 
to non-universal features. At low fields, we can say, at most, 
that the NRG and BA results are consistent within the error- 
bars. This trend is also found in experimental results^ 3 -. In 
experiments where the splitting was found to exceed the pre- 
dicted splitting in the high-field rang o 44,45 , this is likely due 
to the non-universal effects which are expected in systems 
which are not in the extreme Kondo limit, i.e., the regime 
where the Kondo temperature is lower by many orders of the 
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Figure 8: (Color online) The Kondo resonance splitting as a func- 
tion of the external magnetic field. The error bars correspond to the 
"pessimistic" and "optimistic" error estimates. Model parameters are 
U = 0.1, e d = -U/2, T = 0.008 in units of half-bandwidth of the 
(flat) conduction band. The broadening parameter is a — 0.075. The 
data labeled as "Moore Wen" are taken from Fig. 2 in Ref.l37l 



magnitude compared to all other energy scales in the prob- 
lem. Such non-universal effects are expected in the Anderson 
model 3 -^ but also in the Kondo model 4 -^ Furthermore, one 
should take into account the strong asymmetry of the Kondo 
peaks in strong magnetic field: perturbative renormalization 
group calculations 39 show that the maximum of the spectral 
peak is located at cj > g^isB, while the center of the left 
flank of the peak appears to be position almost precisely at 
uj = g/is-B. For very large fields, the peak itself is no longer 
observable and one is left with a step in the spectral function, 
which should be considered as the sole remnant of the Kondo 
resonance. Finally, it should also be noted that for meaning- 
ful comparison with the experimental results, it is necessary 
to take into account the non-equilibrium effects if the splitting 
is extracted from the conductance at finite source-drain volt- 
age in quantum-dot setups with symmetric coupling to both 
leads^. 

We are now in the position to critically discuss the recent 
work on the scaling of the magnetic-field-induced Kondo res- 
onance splitting—, the comment concerning that work^L, and 
the reply offered by the authors of the original work—. In par- 
ticular, it has been claimed 42 that the position of the Kondo 
resonance in the total spectral function does not approach its 
position in the spin-resolved spectra for high magnetic fields, 
in contradiction to what has been found in some previous 
works 4 -^ and that the splitting shows non-universal behavior 
even for modest B /Tjc ratio of order 10. Both conclusions 
have been shown in Ref . [U to be a consequence of the spectral 
function overbroadening due to an excessively large broaden- 
ing parameter a = 0.8 and it was pointed out that different 
results are obtained with smaller broadening a = 0.075. In 
reply, it has been claimed that the value of a = 0.075 is too 
small and that a = 0.4 is a more appropriate choice^. 

The results of the present work, in particular Figs. |6]and[7] 
make it possible to go beyond the purported "certain arbitrari- 
ness" in the choice parameters^ and elucidate to what extend 



the NRG can provide a definitive answer to the problem of 
the magnetic-field-induced Kondo peak splitting. Two points 
need to be emphasized: i) while a = 0.8 is clearly too large (it 
leads to an error in excess of 30% in determining the Kondo 
peak positions) and a = 0.4 is better (15% error), it is crucial 
to go in the a —> limit in order to obtain a result with ac- 
curacy in the percent range, thus a = 0.075 is a good choice 
for N z = 12; ii) using the "pessimistic" error estimate, there 
is a sizable overlap of the confidence regions for a = 0.075 
and a = 0.4, thus there is non-negligible possibility that a cal- 
culation performed for a certain fixed value of N would yield 
similar results for the peak position using both values of a, see 
Fig. |6l although this conclusion is likely to be too pessimistic 
and a different error estimate suggest a clear difference be- 
tween a — 0.4 and a = 0.075 results, see Fig. [7] It may be 
thus concluded that a = 0.075 is a more appropriate choice of 
the broadening parameter and that the results and conclusions 
of Ref. 42 are questionable due to spectral overbroadening. In 
particular, in the a — > limit, the Kondo peak positions in 
the total and spin-resolved spectral functions approach in the 
high-field limi t 31,41 . The conclusion of Ref. [42| that the slope 
coefficient of the ratio of the Kondo peak splitting over mag- 
netic field is 2/3 in the small field limit, as expecte d 34 ' 35 , is 
surprising given the significant shift of the Kondo peak due 
to overbroadening; this result may be simply fortuitous. Us- 
ing log-Gaussian broadening and taking into account the error 
bars in the small-a limit, such slope determination cannot be 
made in a reliable way using the NRG, see Fig. [8] Further- 
more, it may also be remarked that the number of states kept 
in Ref.l42l i.e. N — 150, is too small to obtain well converged 
results (see Fig. |2), irrespective of the broadening procedure 
used. 



VI. INELASTIC (SPIN-FLIP) TUNNELING STEP AND 
THE LARGE-FIELD LIMIT 

If U is much larger than T, for instance U/T = 100, the 
Kondo temperature is for all practical purposes equal to zero, 
since experiments are performed at a finite temperature which 
is, in this case, larger than Tk by orders of magnitude. The 
impurity then behaves much like a free spin, as long as the ex- 
ternal magnetic field is not comparable to the atomic energy 
scales (e, U). In experiments, for example in the inelastic 
spin-flip tunneling spectroscopy using a scanning tunneling 
miscroscope (STM) 4 ^, one can induce inelastic scattering by 
injecting electrons from the STM tip into the adsorbed im- 
purity with energy exceeding the characteristic energy of a 
spin-flip event, i.e., above the Zeeman energy. Due to high 
relevance for STM experiments, it is of substantial interest to 
study the spectral function of the Anderson impurity model 
in the vicinity of the on-set of inelastic scattering. We study 
three aspects of this problem: i) the line-shape of the step in 
the spectral function at the on-set of spin-flip scattering; ii) 
the position of this step as B approaches the atomic scales; 
iii) merging of the step with the atomic peak. 

We express the magnetic field in energy units (Zeeman 
splitting) as b = gusB. For b small compared with the 



7 



<C 0.5 



< 0.5 





1 1 

<C 0.5 


^ 1 

s, 

<C 0.5 



gH R B=10- 



J 



— NRG 

— pert. RG 



gM=10" 



J 



J 



r 



_L 



1 

CO/B 



Figure 9: (Color online) Spectral function of the Anderson impurity 
model in high magnetic field, B 3> Tk ■ Black line corresponds to 
the NRG results, the two ligher (brown onlines) lines demarkate the 
confidence region, while the light (red online) curve corresponds to 
the analytical perturbative RG function from Ref.l39l 



atomic scales, but b ^> Tk, we find that the spectral func- 
tion around uj = b takes the form of a step, see Fig. [9] This 
regime has been studied before using both NRG and pertur- 
bative renormalization group (RG) techniques, see Fig. 5 in 
Ref.[39l Our results are consistent with these studies. Taking 
into account the finite broadening width in the NRG calcula- 
tions, an excellent agreement is found with perturbative RG as 
long as b is much smaller compared to U (see the deviations 
foib/U = 1/5, bottom-most panel in Fig. |9j note that the per- 
turbative RG calculation is performed for the effective Kondo 
model, not for the Anderson model). In experiments, finite 
temperature will play a similar smoothing effect as spectral 
broadening in NRG, thus one can indeed expect to observe a 
step-like spectral line-shape. The step, whose center is always 
located at cj = b, can be interpreted as the on-set of the inelas- 
tic (spin-flip) scattering, which is observed by the spin-flip 
spectroscopy in system that do not exhibit the Kondo effect 4 -. 
It should be noted that there is no discernable peak at u) = b: 
in this b ^> Tk regime, the field-split Kondo resonance has 
become so asymmetric that it takes the form of a relatively 
sharp step^. The width of the step as determined by the NRG 
matches that expected for a unit-step function broadened by 
the kernel, i.e., the intrinsic width of the step is very small, 
presumably equal to T « ■kB/{1& In 2 [B /Tk]), the transverse 
spin relaxation rato^. We also point out that nothing notice- 
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Figure 10: (Color online) Spectral function of the Anderson impurity 
model in very high magnetic field, B 3> Tk and B ~ 17/2. The 
dark lines are the averaged NRG results, while the lighter (brown 
online) lines indicate the confidence region. The discretization errors 
are more pronounced on the high-frequency side of the spectral step 
at w = b = gfisB. 



able happens on the scale b ~ T = 10 3 , nor on the scale 
6- VUT = 10~ 2 . 

We now study the regime where b is comparable to the 
atomic scales, Fig.[l0] We again observe that there is always 
a spectral step exactly at u = B, see right panels in Fig. [9] 
The steps becomes more diffuse as it merges with the atomic 
peak for b > U/2. 

In the atomic limit (V — > 0), the impurity energy levels are 



E 

El 
Eo 



0, 

e + b/2, 
e - 6/2, 
2e + U. 



(8) 



In the particle-hole symmetric point, S = £ + (7/2 = 0, one 
has Ei = Eo, thus for any non-zero value of the magnetic 
field, the state | \) is the ground state and the spin-up spectral 
function has a peak at 

w = E 2 - E x = e + U + b/2 = 5 + U/2 + b/2. (9) 

We find that for b>U, the peak position indeed deviates only 
little from ujq, see left panels in Fig.fTOl 
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VII. CONCLUSION 

An analysis of the spectral functions calculated using the 
NRG technique shows that there is always some variance due 
to the discretization and truncation errors. For the range of 
values of N which is suitable for practical NRG calculations, 
the obtained spectral functions do not converge; instead, the 
variance of the truncation errors appears to be approximately 
constant as a function of TV (for large N). Using a large broad- 
ening parameter does not solve the problem, but merely masks 
it. Furthermore, overbroadening errors appear to be a much 
larger reason for concern than the discretization and trunc a- 
tion errors. Accurate calculations should therefore aim for 
obtaining the results in the limit of zero broadening width, 
taking into account the constraints imposed by the systematic 
NRG errors. Errors should be quantified, not ignored. 

We applied the procedure to estimate the errors in the NRG 
results for the Kondo resonance splitting in the external mag- 
netic field. The systematic errors preclude the study of the 
small-field limit. For intermediate fields, however, it is pos- 
sible to calculate the splitting ratio u>K/g^BB with an esti- 
mated error of a few percent, which is reasonably accurate. 



Good agreement is found with the Bethe Ansatz results for 
the peak splitting by Moore and Wen^ in the regime of small 
and intermediate magnetic fields, where Anderson and Kondo 
models are equivalent. 

To find a definitive quantitative solution of the problem of 
the Kondo resonance splitting in the external magnetic field, it 
will be necessary to devise numerical techniques that can re- 
duce the discretization and truncation artifacts even further. 
Recent developments where the NRG matrix product state 
representation is refined using the density-matrix renormal- 
ization group (DMRG) procedure or variationally^^L may be 
very valuable in this respect since they might allow perform- 
ing calculation at much reduced discretization parameters A 
and optimizing the excited state energies by sweeping. 
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